class BeringEcosystem extends Ecosystem {

  // local variables to store all statevars and the fluxes that go in arrays -------------------------------
  
  Statevar PAR, NH4, NO3, SP, LP, SD, LD, MZ, Eu;
  Statevar[] C, aC; // copepod stages and associated ages
  String[] prey, pred, stages, stage_longnames;
  int Nprey, Npred, Nstages;
  float[][] phi;
  float botNO3_max, botNO3_Hs, surfNO3_max, surfNO3_Hs; // stores a fit to surface and bottom NO3 as a function of water depth
  float small = 1e-6;

  Environment env; // so the model can get access to T, kbot, and dz
  
  
  // model definition ----------------------------------------------------------------------------------------


  BeringEcosystem(String xlsName, String copepodName) {
    addParamsFromXls(this, xlsName, 0); // read in all the free params from an Excel file
    addCopepodFromXls(this, xlsName, copepodName); // reads in more params and the names of the stages
    
    // statevars
    PAR = addGhostVar("PAR","PAR","W/m^2");    
    NH4 = addVar("NH4","NH4","µM");
    NO3 = addVar("NO3","NO3","µM");
    SP = addVar("small phytos","SP","µM");
    LP = addVar("large phytos","LP","µM");
    MZ = addVar("microzooplankton","MZ","µM");
    Eu = addVar("euphausiids","Eu","µM");
    SD = addVar("small detritus","SD","µM");
    LD = addVar("large detritus","LD","µM");
    C = new Statevar[stages.length];
    for (int i=0; i<C.length; i++) {
      C[i] = addVar(stage_longnames[i], stages[i], "µM");
    }
    aC = new Statevar[stages.length-1];
    for (int i=0; i<aC.length; i++) {
      aC[i] = addVar("relative age of " + stages[i], "a"+stages[i], "");
    }

    // phytoplankton growth        
    addFlux("uptake NH4-SP", "Fup_NH4_SP", NH4, SP);
    addFlux("uptake NH4-LP", "Fup_NH4_LP", NH4, LP);
    addFlux("uptake NO3-SP", "Fup_NO3_SP", NO3, SP);
    addFlux("uptake NO3-LP", "Fup_NO3_LP", NO3, LP);
    addFlux("injection of ice algae", "FiceAlgae", null, LP);
    
    // pred-prey lists
    Nstages = stages.length;
    Nprey = 3 + Nstages;
    Npred = 2 + Nstages;
    prey = new String[Nprey];
    prey[0] = "SP";
    prey[1] = "LP";
    prey[2] = "MZ";
    for (int i=3; i<Nprey; i++) prey[i] = stages[i-3];
    pred = new String[Npred];
    pred[0] = "MZ";
    for (int j=1; j<Npred-1; j++) pred[j] = stages[j-1];
    pred[Npred-1] = "Eu";
    
    // grazing preferences
    phi = new float[Npred][Nprey];
    
    // zooplankton grazing
    for (int j=0; j<Npred; j++) {
      for (int i=0; i<Nprey; i++) {
        addFlux("grazing of " + prey[i] + " by " + pred[j], "Fgraz_"+prey[i]+"_"+pred[j], getVar(prey[i]), getVar(pred[j]));
      }
      addFlux("egestion by " + pred[j] + " (" + j + ")",  "Feg_"+pred[j],  getVar(pred[j]), LD); // all predators send egestion to LD, except...
      addFlux("excretion by " + pred[j] + " (" + j + ")",  "Fex_"+pred[j],  getVar(pred[j]), NH4);
    }
    getFlux("Feg_MZ").to = SD; // ... except MZ
    
    // mortality, remineralization, nitrification
    // (sinking is added in the physical environment)
    addFlux("SP mortality", "Fmort_SP", SP, SD);
    addFlux("LP mortality", "Fmort_LP", LP, SD);
    addFlux("MZ mortality", "Fmort_MZ", MZ, SD);
    for (int i=0; i<Nstages; i++) {
      addFlux(stages[i] + " mortality", "Fmort_"+stages[i], C[i], LD);
    }
    addFlux("Eu mortality", "Fmort_Eu", Eu, LD);
    addFlux("SD remineralization", "Fremin_SD", SD, NH4);
    addFlux("LD remineralization", "Fremin_LD", LD, NH4);
    addFlux("nitrification", "Fnitr", NH4, NO3);
    
    // copepod life cycle
    for (int i=0; i<Nstages-1; i++) {
      Flux F = addFlux("molting, " + stages[i] + " to " + stages[i+1], "Fmolt_"+stages[i], C[i], C[i+1]);
      addFlux("aging of " + stages[i], "Faging_"+stages[i], null, aC[i]);
    }
    addFlux(stages[Nstages-1] + " egg production", "FeggProd", C[Nstages-1], C[0]);
    
 
    botNO3_max = getParam("botNO3_max").initial; // if these params are changed on the fly, nothing will happen. If this is changed, note that the i.c. gets set in BeringWaterCol before .current is filled in
    botNO3_Hs = getParam("botNO3_Hs").initial;
    surfNO3_max = getParam("surfNO3_max").initial;
    surfNO3_Hs = getParam("surfNO3_Hs").initial;
  }

  
  void updatePrivateParams() {
    // fill in matrix of prey preferences
    for (int i=0; i<Nprey; i++) {
      for (int j=0; j<Npred; j++) {
        FreeParam phiji = getParam("phi_" + pred[j] + "_" + prey[i]);
        if (phiji != null) { // explicitly set prey preference
          phi[j][i] = phiji.initial;
        } else if (prey[i].equals("LP") && getParam("I0_"+pred[j]).initial > 0) { // everything that eats eats LP; by definition phi_*_LP = 1
          phi[j][i] = 1;
        } else { // j doesn't eat i
          phi[j][i] = 0;
        }
      }
    }
  }
  
  
  // some shortcuts for the flux calculation
  float mu0(String P) {
    return getParam("mu0_" + P).current; 
  }
  float alpha(String P) {
    return getParam("alpha_" + P).current;
  }
  float lightlim(String P, int k) {
//    float LL = 1 - exp(-alpha(P) * PAR.current[k] / mu0(P));
//    float LL0 = 1 - exp(-alpha(P) * getParam("PARcomp").current / mu0(P));
    float LL = alpha(P) * PAR.current[k] / sqrt(sq(mu0(P)) + sq(alpha(P) * PAR.current[k]));
    float LL0 = alpha(P) * getParam("PARcomp").current / sqrt(sq(mu0(P)) + sq(alpha(P) * getParam("PARcomp").current));
    return max(0, LL - LL0) / (1 - LL0);
  }
  float nutlim(String N, String P, int k) {
    float ksmin = getParam("ksmin_" + P + "_NO3").current;
    if (N.equals("NH4")) ksmin *= getParam("rel_ks_NH4").current;
    float Ntot = NH4.current[k] + NO3.current[k];
    return getVar(N).current[k] / (small + ksmin + 2*sqrt(ksmin*Ntot) + Ntot);
  }
  float stageDuration(String stage, int k) {
    float a = getParam("bele_a_" + stage).current;
    float b = getParam("bele_b_" + stage).current;
    float alpha = getParam("bele_alpha_" + stage).current;
    return a * pow(temp(k) + alpha, b);
  }
  float temp(int k) {
    return env.getVar("temp").current[k];
  }
 
 
  // surface and deep nitrate as a function of water-column depth
  // note: if this changes to be a different function of H, use abs(H), or else revisit exactly how H is defined in the forcing files
  float surfaceNO3(float H) {
    return surfNO3_max * sq(H) / (sq(surfNO3_Hs) + sq(H));
  }
  float bottomNO3(float H) {
    return botNO3_max * sq(H) / (sq(botNO3_Hs) + sq(H));
  }
  float NO3bc(float z, float H) { // NO3 b.c., interpolating linearly between the surface and bottom regressions
    float a = constrain(abs(z) / abs(H), 0, 1);
    return a * bottomNO3(H) + (1-a) * surfaceNO3(H);
  }

  void restoreBottomNO3() {
    int kbot = ((Slice)env).kbot[0];
    float H = -((Slice)env).z_rho[kbot];
    float botNO3 = bottomNO3(H);
    getVar("NO3").current[kbot] = botNO3;
    for (int k=0; k<kbot; k++) getVar("NO3").current[k] = 0; // for clarity, set NO3 = 0 below the bottom
  }



  void calcFluxes() {
    // impose a minimum biomass on all living compartments
    float minConc = 0.0001;
    for (int k=0; k<Ncells; k++) {
      SP.current[k] = max(SP.current[k],minConc);
      LP.current[k] = max(LP.current[k],minConc);
      MZ.current[k] = max(MZ.current[k],minConc);
      Eu.current[k] = max(Eu.current[k],minConc);
      for (int i=0; i<Nstages; i++) {
        C[i].current[k] = max(C[i].current[k],minConc);
      }
    }
    
    // food limitation (calculated from max concentration), weighted using prey prefs for each predator
    // the value for MZ is overwritten down below, since MZ are limited by local, not max, food concentration
    float[] foodLim = new float[Npred];
    for (int j=0; j<Npred; j++) {
      foodLim[j] = 0;
      float maxFood = 0;
      for (int k=((Slice)env).kbot[0]; k<Ncells; k++) {
        float foodConc_k = 0;
        for (int i=0; i<Nprey; i++) {
          foodConc_k += phi[j][i] * getVar(prey[i]).current[k];
        }
        if (foodConc_k > maxFood) maxFood = foodConc_k;
      }
      foodLim[j] = maxFood / (small + getParam("K0_" + pred[j]).current + maxFood);
    }
    
    restoreBottomNO3();

    getFlux("FiceAlgae").current[Ncells-1] = - getParam("IAint").current * min(0, ((BeringWaterCol)env).ddt_icecover()) / ((Slice)env).dz[Ncells-1];
    
    for (int k=((Slice)env).kbot[0]; k<Ncells; k++) {
      // uptake
      float TFactor = pow(getParamVal("Q10_P"), temp(k) / 10.0);
      getFlux("Fup_NH4_SP").current[k] = nutlim("NH4","SP",k) * lightlim("SP",k) * mu0("SP") * TFactor * SP.current[k];
      getFlux("Fup_NO3_SP").current[k] = nutlim("NO3","SP",k) * lightlim("SP",k) * mu0("SP") * TFactor * SP.current[k];
      getFlux("Fup_NH4_LP").current[k] = nutlim("NH4","LP",k) * lightlim("LP",k) * mu0("LP") * TFactor * LP.current[k];
      getFlux("Fup_NO3_LP").current[k] = nutlim("NO3","LP",k) * lightlim("LP",k) * mu0("LP") * TFactor * LP.current[k];

      // mortality & regeneration
      getFlux("Fmort_SP").current[k] = getParam("m_SP").current * SP.current[k];
      getFlux("Fmort_LP").current[k] = getParam("m_LP").current * LP.current[k];
      getFlux("Fmort_MZ").current[k] = getParam("m_MZ").current * sq(MZ.current[k]);
      for (int i=0; i<stages.length; i++) {
        getFlux("Fmort_"+stages[i]).current[k] = getParam("m_"+stages[i]).current * C[i].current[k];
      }
      getFlux("Fmort_Eu").current[k] = getParam("m_Eu").current * Eu.current[k];
      TFactor = pow(getParamVal("Q10_regen"), temp(k) / 10);
      getFlux("Fremin_SD").current[k] = getParam("R_remin").current * TFactor * SD.current[k];
      getFlux("Fremin_LD").current[k] = getParam("R_remin").current * TFactor * LD.current[k];
      getFlux("Fnitr").current[k] = getParam("R_nitr").current * NH4.current[k];
      
      // grazing
      float TFactorMeso = pow(getParamVal("Q10_meso"), temp(k) / 10.0);
      for (int j=0; j<Npred; j++) {
        getFlux("Feg_" + pred[j]).current[k] = 0;
        getFlux("Fex_" + pred[j]).current[k] = 0;
        float epsilon, fNH4;
        float foodConc_jk = 0;
        for (int i=0; i<Nprey; i++) {
          foodConc_jk += phi[j][i] * getVar(prey[i]).current[k];
        }
        if (pred[j].equals("MZ")) {
          epsilon = getParam("epsilon_MZ").current;
          fNH4 = getParam("f_MZ_NH4").current;
          foodLim[j] = foodConc_jk / (small + getParam("K0_MZ").current + foodConc_jk);
          TFactor = pow(getParamVal("Q10_MZ"), temp(k) / 10.0);
        } else {
          epsilon = getParam("epsilon_meso").current;
          fNH4 = getParam("f_meso_NH4").current;
          TFactor = TFactorMeso;
        }
        for (int i=0; i<Nprey; i++) {
          float preyFrac_ijk = phi[j][i] * getVar(prey[i]).current[k] / (small + foodConc_jk);
          float Fgraz_ijk = getParam("I0_"+pred[j]).current * getVar(pred[j]).current[k] * preyFrac_ijk * foodLim[j];
          Fgraz_ijk *= TFactor;
          getFlux("Fgraz_"+prey[i]+"_"+pred[j]).current[k] = Fgraz_ijk;
          getFlux("Feg_"+pred[j]).current[k] += (1 - epsilon - fNH4) * Fgraz_ijk;
          getFlux("Fex_"+pred[j]).current[k] += fNH4 * Fgraz_ijk;
        }
      }
      
      // copepod life cycle
      float delta = 0.1;
      float u0;
      for (int i=0; i<Nstages-1; i++) {
        u0 = 1. / stageDuration(stages[i],k);
        getFlux("Fmolt_"+stages[i]).current[k] = u0 / delta * (0.5 + 0.5*tanh((aC[i].current[k] - 1) / delta)) * C[i].current[k];
        getFlux("Faging_"+stages[i]).current[k] = u0 - getFlux("Fmolt_"+stages[i]).current[k] * aC[i].current[k] / (small + C[i].current[k]);
      }
      float netgrowth = getFlux("Fex_"+stages[Nstages-1]).current[k] / getParam("f_meso_NH4").current * getParam("epsilon_meso").current; // inelegant formula for net adult growth. Easier than summing over all the grazing fluxes
      getFlux("FeggProd").current[k] = netgrowth; // all adult net growth goes into egg production
    }       
  } 
}


